# const programs2012 = let
#     f = matopen(datapath*"/"*"additionalProgramInfo_2012.mat")
#     programdata = MAT.read(f,"Programs2012")
#     close(f)
#     programdata
# end
function mymean(x)
    xtrim = x[x.>0]
    isempty(xtrim) ? 0 : mean(xtrim)
end

#get program bins w/ equal numbers of students
const programselectivitybin = let
    (cc,lv,temps) = dataset[2012]
    (cf_constants,_) = counterfactuals[2012]
    J = length(cc.programIDs)
    selectivity = [mymean( cf_constants.priorityEligible[(cc.enroll.==jj) .& .!(cc.specialAdmit), jj]) for jj=1:J]
    selectivity_rank = sortperm(vec(selectivity))
    totalEnrollment = [sum(cc.enroll .== j) for j=1:J]
    enr = cumsum(totalEnrollment[selectivity_rank])
    cuts = maximum(enr) .* (.1:.1:.9)
    mygroup = ones(Int,length(selectivity_rank))
    for cutpoint in cuts
        mygroup[enr .> cutpoint] .+= 1
    end
    programdecile = mygroup[invperm(selectivity_rank)]
end

#fill in lv
let
    t = save_lv[end]
    year=2012
    JLD2.@load(outputpath*"/latentvariables$(year)_$t.mat",lv_out)
    (cc,lv,temps) = dataset[2012]
    lv.U .= lv_out.U
    lv.U0 .= lv_out.U0
    lv.rc .= lv_out.rc
    lv.Astar .= lv_out.Astar
    lv.H .= lv_out.H
end

#run counterfactuals
c_means = Any[]
c_individuals = Any[]
c_byIter = Any[]
for (r,t) in enumerate(save_lv)
    year=2012
    println("running counterfactuals: iter $t")
    JLD2.@load(outputpath*"/latentvariables$(year)_$t.mat",lv_out)
    theta_t = history[findfirst(save_params.==t)]
    (df_means, df_individuals, df_byIter) = runCounterfactuals(theta_t,dataset,counterfactuals,lv_out)
    push!(c_means,df_means)
    df_individuals[!,:draw] = r .* ones(Int,nrow(df_individuals))
    push!(c_individuals,df_individuals)
    push!(c_byIter,df_byIter)
end
df = vcat(c_individuals...)

# c_individuals_extra = Any[]
# for (r,t) in enumerate(save_lv)
#     year=2012
#     println("running more counterfactuals: iter $t")
#     JLD2.@load(outputpath*"/latentvariables$(year)_$t.mat",lv_out)
#     theta_t = history[findfirst(save_params.==t)]
#     (df_means, df_individuals, df_byIter) = runExtraCounterfactuals(theta_t,dataset,counterfactuals,lv_out)
#     df_individuals[!,:draw] = r .* ones(Int,nrow(df_individuals))
#     push!(c_individuals_extra,df_individuals)
# end
# df_extra = vcat(c_individuals_extra)

# CSV.write(outputpath*"/counterfactuals_long.csv",df)
# varnames = [join((split(x,"_")[1:end-1]),"_") for x in filter(
#     contains("enroll"),names(df))]
# for name in varnames
#     (cf_constants,_) = counterfactuals[2012]
#     n_enroll = Symbol(name*"_enroll")
#     df[!,Symbol("$(n_enroll)_Any")] = df[!,n_enroll].>0
#     df[!,Symbol("$(n_enroll)_G25")] = [e > 0 && !(cf_constants.G8[e]) for e in df[!,n_enroll]]
#     df[!,Symbol("$(n_enroll)_Inefficient")] = df[!,n_enroll] .!= df[!,:zeroalpha_enroll]
#     #
#     for key in ["welfare","graduate"]
#         df[!,Symbol(name*"_"*key*"_Any")] = df[!,Symbol(name*"_"*key)] .* df[!,Symbol("$(n_enroll)_Any")]
#         df[!,Symbol(name*"_"*key*"_G25")] = df[!,Symbol(name*"_"*key)] .* df[!,Symbol("$(n_enroll)_G25")]
#         df[!,Symbol(name*"_"*key*"_Inefficient")] = df[!,Symbol(name*"_"*key)] .* df[!,Symbol("$(n_enroll)_Inefficient")]
#     end
# end

# keepcols = getindex(names(df), [occursin("Any",x) || occursin("G25",x) || occursin("Inefficient",x) for x in names(df)])
df_small = combine(groupby(df,:ind),[x=>mean=>x for x in names(df)[3:end]])
df_small[!,:type] = dataset[2012][1].typ
#
# for col in names(df_small)
#     if eltype(df_small[!,col]) == Float64
#         df_small[!,col] = round.(df_small[!,col],digits=5)
#     end
# end

for col in names(df)
    if eltype(df[!,col]) == Float64
        df[!,col] = round.(df[!,col],digits=4)
        df_small[!,col] = round.(df_small[!,col],digits=4)
    end
end
CSV.write(outputpath*"/counterfactuals_big.csv",df)
CSV.write(outputpath*"/counterfactuals_means.csv",df_small)
df_bytype = vcat( [(df = copy(cm); df[!,:draw] = [t,t,t,t,t]; df) for (t,cm) in enumerate(c_means)]... )
CSV.write(outputpath*"/counterfactuals_bytype.csv",df_bytype)

#JLD2.@save(datapath*"/counterfactualResults.jld2",c_means,c_individuals,c_byIter)

#=
#avg. welfare of people enrolling in G25 under "alpha=0" counterfactual
sum( df_small[!,:zeroalpha_welfare_G25])/sum(df_small[!,:zeroalpha_enroll_G25])
#avg. welfare of people enrolling in G25 under "as if 2011" counterfactual
sum( df_small[!,:asif2011_welfare_G25])/sum(df_small[!,:asif2011_enroll_G25])
#avg welfare loss, relative to efficient, among people who got a different allocation
sum( df_small[!,:zeroalpha_welfare_Any] .- df_small[!,:asif2011_welfare_Any]) /
    sum(df_small[!,:asif2011_enroll_Inefficient]) #note numerator is zero for ppl whose allocation didnt change
sum(df_small[!,:asif2011_enroll_Inefficient]) #number of such people
=#
